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ABSTRACT 

The highly pathogenic influenza strains H5N1 and 
H1N1 are currently treated with inhibitors of the viral 
surface protein neuraminidase (Nl). Crystal structures 
of Nl indicate a conserved, high affinity calcium bind- 
ing site located near the active site. The specific role of 
this calcium in the enzyme mechanism is unknown, 
though it has been shown to be important for enzy- 
matic activity and thermostability. We report molecular 
dynamics (MD) simulations of calcium-bound and cal- 
cium-free Nl complexes with the inhibitor oseltamivir 
(marketed as the drug Tamiflu), independently using 
both the AMBER FF99SB and GROMOS96 force fields, 
to give structural insight into calcium stabilization of 
key framework residues. Y347, which demonstrates 
similar sampling patterns in the simulations of both 
force fields, is implicated as an important Nl residue 
that can "clamp" the ligand into a favorable binding 
pose. Free energy perturbation and thermodynamic 
integration calculations, using two different force 
fields, support the importance of Y347 and indicate a 
+3 to +5 kcal/mol change in the binding free energy 
of oseltamivir in the absence of calcium. With the im- 
portant role of structure-based drug design for neur- 
aminidase inhibitors and the growing literature on 
emerging strains and subtypes, inclusion of this cal- 
cium for active site stability is particularly crucial for 
computational efforts such as homology modeling, vir- 
tual screening, and free energy methods. 
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INTRODUCTION 

Type A influenza virus is becoming a worldwide pandemic 
threat due to its virulence and transmissibility in people. The 
viral surface protein neuraminidase, along with hemagglutinin, 
classifies influenza subtypes and fulfills an important role in vi- 
ral propagation by cleaving a terminal sialic acid from host cell 
surfaces; thus, neuraminidase is a key drug target for control- 
ling flu infection. The enzyme has two distinct families, ^ group 
1 (Nl, N4, N5 and N8) and group 2 (N2, N3, N6 and N7), 
with group 1 including the Nl subtype found in recent out- 
breaks of the highly pathogenic subtype H5N1, as well as the 
1918 epidemic and "swine flu" H1N1. The protein exhibits a 
high rate of antigenic shift and drift, yielding proteins with 
heterogeneous sequences that can reduce sensitivity to the clin- 
ically used inhibitors oseltamivir (marketed as Tamiflu by 
Roche Pharmaceuticals, Basel, Switzerland), zanamivir (mar- 
keted as Relenza by GlaxoSmithKline, Bentford), and peramivir 
(BioCryst Pharmaceuticals, Birmingham, AL). 

A high-affinity calcium binding site near the active site has 
been identified through crystallography^ - ^ and Proton 
Induced X-ray Emission (PIXE) experiments,^ though the cal- 
cium density was not resolved or discussed for some holo 
structures of group 1 neuraminidase. ^ The specific role of this 
calcium ion is unknown, but experiments have supported 
its importance for wild-type enzyme activity^ and enzyme 
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thermostability'. Commonly performed fluorometric 
neuraminidase activity assays include calcium salts in the 
activity buffer and follow the protocol developed origi- 
nally by Potier et al.fi where increased calcium ion con- 
centration was found to augment activity. In addition, a 
recent crystal structure of calcium-deficient group 2 
neuraminidase shows destabilization of the key active site 
residue R292 near the calcium binding site.^ 

Although dynamics of Nl ligand binding have been 
examined computationally,^ - ^ no dynamic, structural 
study has investigated the contribution of this calcium to 
neuraminidase substrate binding. Here, we compare sim- 
ulations of calcium-bound and calcium-free complexes of 
Nl neuraminidase bound to the FDA-approved, nano- 
molar- affinity drug oseltamivir. To quantify the impact of 
this calcium on oseltamivir binding affinity, we calculate 
the free energy of binding for both complexes. This study 
also checks for force field dependency of the observed 
calcium effects on Nl dynamics through comparison of 
simulations and free energy calculations performed with 
the AMBER FF99SB and GROMOS96 force fields. 

We emphasize the impact of calcium binding on integ- 
rity of the active site and accurate calculation of binding 
free energies. In particular, for Nl, calcium stabilizes 
Y347 for interaction with the ligand and helps to main- 
tain an optimal binding pose with other Nl binding site 
residues. Incorporation of calcium as an integral part of 
the active site is particularly crucial given the importance 
of structure-based drug design for Nl inhibitor develop- 
ment to date^ and prevalance of homology models^ 5 '^ 
for emerging subtypes and resistance mutations, virtual 
screening, 1'>1° and free-energy methods 19-21 f or m ; s 
protein. 



MATERIALS AND METHODS 

Nl monomer simulations were performed using the 
GROMOS05 software for biomolecular simulation 22 and 
the GROMOS96 force field (45A3 parameter set). 23 Mo- 
lecular Dynamics (MD) set up is described in detail in 
the Supporting Information Methods section. Both com- 
plexes of Nl -oseltamivir, with and without the bound 
calcium ion, were simulated in 10 independent trajecto- 
ries, each 4 ns, for a total of 40 ns of simulation for each 
complex. Comprising the 10 simulations were five simu- 
lations generated from the chain B monomer in the osel- 
tamivir-bound Loop 150 "open" crystal structure 
(PDBID: 2HU0) and five simulations started from the 
chain A monomer in the holo Loop 150 "closed" crystal 
structure (PDBID: 2HU4). As the calcium density was 
not present in these crystal structures, overlap with the 
apo Nl structure 2HTY aided in positioning of the ion 
in the protein; the calcium was parametrized in the clas- 
sical force field. To generate the independent trajectories, 
these structures were each initialized with random veloc- 



ities assigned from a Maxwell-Boltzmann distribution at 
5 K. 

Two 100 ns Nl tetramer simulations with the AMBER 
FF99SB force field 24 ' 25 were performed, each with 
atomic coordinates taken from the holo, open Loop 150 
crystal structure (2HU0) and with the calcium inserted 
from overlap with the apo 2HTY structure. Calcium- 
bound simulations used the PMEMD module in AMBER 
lO, 2 ^ whereas calcium-free simulations were performed 
using the Desmond Molecular Dynamics package devel- 
oped by D.E.Shaw Research. 2 '' For details of these simu- 
lations, please see the Supporting Information Methods 
section. Although there are a few differences in the spe- 
cifics of the MD engines used to run the AMBER 
FF99SB simulations (Supporting Information Methods), 
the similarity in RMSF (Supporting Information Fig. 1) 
indicates similar sampling of conformational space. 
AMBER FF99SB trajectories for each monomer of the 
tetramer were extracted and concatenated to approximate 
400 ns of monomer Nl sampling. 

Analysis tools in Visual Molecular Dynamics (VMD) 2 ^ 
were applied for calculation of root mean squared fluctu- 
ation (RMSF) and deviation (RMSD), and monitoring of 
torsion, hydrogen bonds, and salt bridge distances. 
Hydrogen bond criteria applied a donor-acceptor dis- 
tance cutoff of 3.5 A and a 120° cutoff for the donor- 
hydrogen-acceptor angle. For analysis, structures were 
extracted every 10 ps for each 400 ns AMBER FF99SB 
trajectory and 2 ps for the 40 ns GROMOS96 trajecto- 
ries. These structures were used for the GROMOS RMSD 
clustering algorithm, 2 ^ applied using the cluster2 pro- 
gram in GROMOS05 software. 22 A rotational and trans- 
lational fitting was applied to all the C a carbons, fol- 
lowed by application of the RMSD clustering to all atoms 
of 41 residues that comprise the binding site (residues 
117-119, 133-138, 156, 178-179, 196-200, 223-228, 243- 
247, 276, 277, 292-295, 344-347, 371, 401, 402). This 
selection is similar to that used for relaxed complex 
docking of Nl in Cheng et al., but here we exclude 
residues in the highly flexible Loop 150 (residues 147- 
152) and Loop 430 (430-439) regions to focus the algo- 
rithm on the active site portion near the calcium binding 
site. After screening a variety of RMSD cutoff values for 
cluster generation, using the total number and diversity 
of clusters as criteria, a value of 1.5 A was selected for 
the GROMOS96 trajectories, and a value of 1.2 A for 
AMBER FF99SB trajectories. The discrepancy in the cho- 
sen RMSD cutoff, as well as the varied cluster popula- 
tions, can be attributed to the varied potential energy 
landscapes produced by the different force fields, as well 
as the use of multiple, shorter GROMOS96 trajectories 
and fewer, longer AMBER FF99SB simulations. 

Two different, rigorous binding free energy calculations 
were used: free energy perturbation (FEP), with the 
AMBER FF99SB force field, and thermodynamic integra- 
tion (TI) with the GROMOS96 force field. Both bound 



2524 PROTEINS 



N1 Influenza Neuraminidase Dynamics and Binding Free Energy 



and unbound calculations were performed, to give a 
AGp rote i n from decoupling the bound oseltamivir from 
Nl, and a AG water from decoupling of unbound oseltami- 
vir, free in a box of water. Following the thermodynamic 
cycle involved in double decoupling,^'^ 1 AG prote i n is 
subtracted from AG wate r to give AGbind- 

The GROMOS05 22 software was used for both oselta- 
mivir-bound and oseltamivir-unbound TI calculations, 
repeated in independent trajectories (IT-TI). For the 
bound state decoupling, six calculations decoupled osel- 
tamivir from the ion-bound complex, and six decoupled 
the ligand from the ion-free complex, for a total of 12 
Nl TI calculations. For each of these complexes, the cal- 
culations were performed in triplicate, with different 
starting structures that capture diverse conformations of 
Loop 150; three calculations were initialized from a 
closed Loop 150 Nl conformation, and three were initial- 
ized from an open Loop 150 conformation. These start- 
ing structures were extracted from the independent MD 
trajectories described above, after 1 ns of equilibration. 
Six independent oseltamivir-unbound TI calculations 
decoupled oseltamivir from a 37 X 37 X 37 A 3 box of 
SPC water. For each TI, 26 windows span the X range 
from 0.00 (ligand coupled) to 1.00 (ligand decoupled), in 
0.04 intervals. For ligand-protein decoupling, a 500 ps 
equilibration period preceded a 500 ps sampling period, 
from which W- data was recorded; the unbound decou- 
pling calculations involved 400 ps for both equilibration 
and sampling periods per X window. 

To maintain a defined bound state in the protein cal- 
culations, a 0.6 kcal/mol/A~ harmonic restraint constant 
was applied to the C4 atom of oseltamivir and the asso- 
ciated volume correction 3 * + RTln(CV) for the free 
energy was calculated to be —5.5 kcal/mol. Soft-core 
interaction potentials 32 (5 LJ = 0.5 and s c = 0.5) were 
used for ligand atoms involved in the perturbation to 
avoid end point singularities and to enhance sampling. 
For each X window, the statistical inefficiency of the raw 
W data was calculated using the method described and 
coded by Chodera et a\?^ The decorrelated data was 
averaged to approximate (IP), which was integrated by 
the trapezoidal rule over the 26 X windows for each of 
N = 6 independent calculations, as: 



Subtraction of the AG prot em from AG wate r gives a 
AGbind- A standard error -S= was computed for both 
AGp rote i n an d AG water over the N = 6 independent runs 
and propagated for AGb; n d (listed in Table III). 

To examine the advantage of shorter, independent tra- 
jectories over a single, long trajectory in the TI calcula- 
tions, one run for each set of protein TI calculations was 
selected for the extension of equilibration and sampling 



times (1 ns equilibration and 3.8 ns sampling per X), for 
an 85.8 ns increase in total sampling time for the calcula- 
tion. The single AG prote i n produced from each extended 
run differed from that of the "short" calculation result 
by <4%; however, the associated AGb m d was +3.3 kcal/ 
mol less favorable than the experimentally derived target 
value. This discrepancy is larger than that for the IT-TI 
result for AGbi n d> which is off the target value by +1.3 
kcal/mol (Table III). The extended ion-free result was 
also significantly more unfavorable than the IT-TI value. 
This indicates that these single, extended simulations 
were unable to overcome sampling barriers to generate a 
complete ensemble average {W). Combining multiple TI 
calculations initialized from diverse conformations 
improves sampling and can result in a more correct 
Ar 12,34,35 

^Gjind- 

Free energy calculations on the tetrameric complex 
were performed with the AMBER FF99SB force field in 
Desmond using FER For both the decoupling of 
unbound oseltamivir (to calculate AG water ) and bound 
oseltamivir (for AG prote i n ), 21 windows were used, one 
with full non-bonded interaction, 10 for the annihilation 
of columbic interactions (with full van der Waals interac- 
tions), and 10 with the removal of van der Waals interac- 
tions (with no electrostatics). A softcore potential was 
used with a = 0.5. Decoupling of unbound oseltamivir 
in water was performed on a system of the ligand in a 24 
X 28 X 29 A 3 TIP3P water box with 250 ps of equilibra- 
tion followed by 1 ns of sampling in each window. 
Bound oseltamivir decoupling calculations were per- 
formed on the tetrameric complex (from 2HU0, 
described above) with all four ligand molecules simulta- 
neously removed to enhance sampling. Following 5 ns of 
equilibration of the starting structure, each X window 
was individually equilibrated for 1 ns and then sampled 
for 2 ns. 

Three sets of protein-ligand calculations were per- 
formed, two with the protein coordinates from the 2HU0 
crystal structure (one with the bound calcium and one 
without) and another without calcium and initiated from 
clustering results in which Y347 is flipped out of the 
pocket (see below). Reported work functions were post- 
processed with the Multistate Bennett Acceptance Ratio 
(MBAR). 36 Ener gies were decorrelated based on their 
statistical inefficiencies, as in the GROMOS96 TI simula- 
tions, and error bars representing the analytic errors 
developed in the MBAR implementation, here repre- 
sented by a, were calculated with code made publicly 
available by Shirts and Chodera.^ The free energy of 
decoupling the four ligands from the Nl tetramer was 
divided by four to obtain an average binding free energy 
to the Nl monomer, with a standard error -5- calculated 

V N 

using N = 4. 
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Nl Residue Index 



Figure 1 

Ca RMSF differences for ion-free and ion-bound simulations. 
GROMOS96 (gray) and AMBER FF99SB (black) ion-bound RMSF is 
subtracted from the ion-free simulation RMSF. Asterisks indicate 
residues that contact the calcium ion in the 2HU0 structure. 

RESULTS AND DISCUSSION 

RMSF 

The structure of the Nl calcium binding site includes 
close ion contacts with the backbone carbonyls of D293, 
N294, G297, G345, A346, Y347, and the carboxyi groups 
of both D293 and D324, with similar interactions 
observed in group 1 and 2 crystal structures. i- 4>9 There- 
fore, removal of the ion can be expected to destabilize 
these residues and the loops on which they are located, 
leading to increased fluctuation. Both GROMOS96 and 
AMBER FF99SB ion-bound MD simulations were 
checked for similar RMSF per residue C a (Supporting In- 
formation Fig. 1), with the most flexibility observed for 



Loop 150, Loop 430, and portions of a long disordered 
loop (residues 327-348) surrounding the calcium binding 
site. Plots of the difference in ion-free and ion-bound 
simulation RMSF (Fig. 1) reveal changes in backbone 
RMSF for these regions. Both AMBER FF99SB and GRO- 
MOS96 simulations indicate increased fluctuation near 
the calcium-contacting residues (see peaks near asterisks 
in Fig. 1), with some of the peaks corresponding to resi- 
dues very close to the active site. RMSF changes are 
color-mapped onto the Nl monomer structure in Sup- 
porting Information Figure 2, for an overall view of 
where these changes occur in the protein. Loop 150 and 
430 (labeled in Fig. 1) have changes in flexibility that are 
inconsistent in comparison of the GROMOS96 and 
AMBER FF99SB simulations. 

Y347 torsion 

Of all residues with increased mobility in the Nl ion- 
free simulations, Y347 has the most direct impact on 
ligand binding. The Y347 backbone carbonyl coordinates 
calcium in the Nl crystal structures ^' [pink conforma- 
tion in Fig. 2(a)], stabilizing the residue for interaction 
with the conserved ligand carboxyi group. In Figure 2(b), 
distributions for the tyrosine %l torsion indicate a 
change of Y347 conformation in the ion-free simulations 
(dashed vs. solid lines). Torsions within the range 30- 
105° designate conformations that direct the sidechain 
into the active site for hydrogen bonding to the ligand 
("in" conformation), whereas other torsion ranges corre- 
spond to conformations flipped out of the pocket and 
inaccessible to the ligand ("out" conformation). With the 




Figure 2 

Y347 Sampling in the Nl Active Site, (a) The Nl monomer in complex with oseltamivir (gray) and with bound calcium (gray), as well as the Y347 
"in" conformation (pink) and two representative "out" conformations (yellow and cyan), (b) Normalized probability distributions for the Y347 %\ 
sidechain torsion plotted for ion-bound AMBER FF99SB (black, solid lines) and GROMOS96 (red, solid lines) simulations, compared to the 
corresponding distributions for ion-free simulations (dashed lines). 
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Table I 

Change in % Occurrence of Nl-Oseltamivir Hydrogen Bonds in Ion-Free Simulations Relative to Ion-Bound Simulations 



Carboxyl Ammonium Acetamide 

Ligand moiety R118 R292 R371 Y347 El 19 E227 E277 E227 E277 

AMBER FF99SB GR0M0S96 +6 -13 -13 -14 -30 None a None None None 

-8 -40 -21 -27 -18 -49 -35 -16 -3 



indicates no hydrogen bond observed to meet criteria described in Methods. 

loss of its stabilizing backbone contact with calcium, 
Y347 increases sampling of "out" conformations [see 
cyan and yellow conformations in Fig. 2(a)]. 

Integration of the curves in Figure 2(b) allows for the 
estimation of the total time spent in the "out" Y347 con- 
formations, whereas tracking the torsion time series indi- 
cates how often transitions occur (Supporting Informa- 
tion Fig. 3). In the ion-free simulations, the destabilized 
Y347 flips out of the pocket for 54% of the GROMOS96 
simulation, and for 29% of the AMBER FF99SB simula- 
tion. Examination of each 100-ns monomer AMBER 
FF99SB trajectory reveals an initial "flip" after 5 ns in the 
chain B simulation, but no flip is seen before 13, 28, and 
58 ns for chains A, C, and D, respectively (see Support- 
ing Information Fig. 3). Five of the 10, 4 ns GROMOS96 
ion-free simulations were initialized in a Y347 "out" con- 
formation, four of which do not visit the "in" conforma- 
tion. Those initialized with the "in" Y347 conformation 
flipped within 1.4 ns, whereas one run remained 
unflipped throughout 4 ns (Supporting Information Fig. 
3). Thus, although the varied AMBER FF99SB and GRO- 
MOS96 simulation lengths prohibit a direct comparison 
of duration for the Y347 flipped conformation, the GRO- 
MOS96 simulations seem to allow more facile Y347 flip- 
ping than do the AMBER FF99SB simulations. 

The "in" ligand-accessible conformation is dominant 
in the ion-bound simulations, with the "out" conforma- 
tion occurring for only 20% and 4% of the GROMOS96 
and AMBER FF99SB simulations, respectively. We note 
that not all conformations in the "in" torsion range 
result in hydrogen bond formation, but the Y347 phenol 
is directed toward the ligand carboxyl for potential inter- 
action. 



Hydrogen bonds 

The increased sampling of Y347 "out" conformations 
in the ion- free simulations reduces the availability of this 
residue for hydrogen bonding to oseltamivir. Hydrogen 
bond analysis using geometric criteria (see Methods) 
reveals an overall reduction in the occurrence of hydro- 
gen bonds formed between the protein and ligand for 
the ion-free simulations (Table I). Y347 hydrogen bond 
occurrence is reduced by 27% in the GROMOS96 ion- 
free simulations and by 14% in the AMBER FF99SB ion- 
free simulations. 



Calcium removal affects the stability of other key bind- 
ing site residues as well, resulting in reduced occurrence 
of hydrogen bonds to the ligand. The arginine triad, an 
important neuraminidase binding motif that includes 
R118, R292, and R371, interacts with the carboxyl of 
ligands. ^ A large loss of ligand hydrogen bonds to R292 
and R371 is observed in the ion-free GROMOS96 simu- 
lations, and less severe but significant reductions are also 
observed for these residues in the AMBER FF99SB ion- 
free simulations. In the AMBER FF99SB simulations of 
both complexes, R118 only weakly meets hydrogen bond- 
ing criteria for the ligand. 

With further comparison of changes in hydrogen bond 
occurrence, a difference in sampling becomes evident for 
the simulations with two different force fields. A signifi- 
cant reduction in hydrogen bond occurrence between the 
oseltamivir ammonium and the El 19 carboxyl group is 
observed for both force fields, but the GROMOS96 simu- 
lations favor interactions of the ligand ammonium with 
E227 and E277; in the AMBER FF99SB simulations this 
ligand moiety is consistently oriented toward El 19. 
Hydrogen bonds to residues on the flexible Loop 150 are 
also sampled differently in the AMBER FF99SB and 
GROMOS96 simulations (Supporting Information Table 
I). This inconsistency is also seen in the RMSF changes 
in Figure 1, and could be due to the different MD proto- 
cols used for GROMOS96 and AMBER FF99SB simula- 
tions, as each monomer of the AMBER FF99SB simula- 
tion was initialized in the same open Loop 150 confor- 
mation, but allowed to run for much longer periods 
compared to the GROMOS96 simulations (see Methods). 
In the AMBER FF99SB ion-free simulation, ligand hydro- 
gen bonds to Loop 150 residues D151 and R152 are sig- 
nificantly reduced (Supporting Information Table I) and 
Loop 150 is observed to sample primarily open states 
(Supporting Information Fig. 4). 

The protein-ligand interactions described in Table I 
are also largely electrostatic; monitoring changes in the 
distance between the center of mass for these reacting 
moieties [Supporting Information Fig. 5(a-d)] confirm 
the breaking of salt bridges and corroborate the hydrogen 
bonding data. 

Clustering and structural changes 

RMSD clustering reduces the simulation structural 
data to a set of representative conformations for analysis. 
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Table II 

Summary of Cluster Populations 



GR0M0S9B AMBER FF99SB 



Ion No Ion Ion No Ion 



% % % % 



Cluster 8 


Ensemble 


Y347 


Ensemble 


Y347 


Ensemble 


Y347 


Ensemble 


Y347 


1 


24 


in 


39 


out 


93 


in 


50 


in 


2 


12 


in 


11 


in 






18 


out 


3 


9 


in 


10 


in 






15 


out 


4 


8 


out 


9 


out 






13 


out 


5 


8 


in 


9 


in 










6 


8 


in 


9 


out 










7 


8 


out 


8 


out 










8 


7 


out 















^Clusters listed combine to represent at least 85% of the total ensemble of configurations from the simulation. 



Clusters that together capture at least 85% of the total 
simulation ensemble are listed in Table II, with the con- 
formation of Y347 indicated. The varied cluster popula- 
tions in Table II also highlight differences in the AMBER 
FF99SB and GROMOS96 simulations. The cutoffs were 
chosen to balance resolution of unique clusters with 
reduction of the simulation data to a manageable num- 
ber of representative conformations. The higher selected 
cutoff for the GROMOS96 simulations as well as the 
increased total number of clusters show that both GRO- 
MOS96 simulations allow more fluctuations and 
increased heterogeneity in sampling of conformational 
space; this was also indicated by higher dimensionality 
for GROMOS96 essential dynamics found in a broad 
force field comparison study. ^ The difference in sam- 
pling could also be attributed to the fact that the GRO- 
MOS96 simulations, which consist of ten 4 ns trajecto- 
ries, may benefit from increased sampling due to the use 
of multiple, shorter trajectories,^'^ as opposed to the 
AMBER FF99SB simulations, which consist of fewer 
(four), longer (100 ns) trajectories. 

Structural differences observed by comparison of clus- 
ter central members for the different force fields largely 
correspond to sampling differences observed in hydrogen 
bonding data (Table I) and RMSF changes (Fig. 1). How- 
ever, both AMBER FF99SB and GROMOS96 ion-free 
clusters indicate changes in position of Y347 and N294, 
which both have backbone carbonyls directed toward the 
calcium and have sidechains that extend into the active 
site [Fig. 3(a,c)]. Other common differences occur on a 
solvent exposed loop near the rim of the binding site, in 
the 247-249 region (Fig. 1 and Supporting Information 
Fig. 2). Figure 3 has examples of conformational changes 
represented in the ion-free GROMOS96 and AMBER 
FF99SB clusters, which could be most pertinent for accu- 
racy of docking or other computational free energy 
methods. 

Destabilization of N294 due to its lost backbone con- 
tact to the calcium enables the N294 polar sidechain to 



shift closer into the hydrophobic subpocket^O that 
accepts the oseltamivir aliphatic tail [Fig. 3(a,b)]. This 
shift is further revealed by tracking the proximity of the 
N294 sidechain to the oseltamivir aliphatic tail through- 
out all simulations [Supporting Information Fig. 5(e,f)]. 
Instability of this residue, as well as Y347, has been pro- 
posed as the source of resistance to oseltamivir conferred 
by the N294S mutation in Nl. 4 In the GROMOS96 sim- 
ulations, changes in R292 and E276 conformations also 
bring these highly polar residues closer to the hydrohobic 
tail of oseltamivir [Supporting Information Fig. 5(e)]. 

Another significant conformational change is observed 
in cluster four for the GROMOS96 ion-free simulation. 
Here, Y347 is flipped out of the pocket, and the nearby 
residue R292 retreats away from its electrostatic interac- 
tion with the oseltamivir carboxyl to interact with flank- 
ing glutamates E276 and E277 [Fig. 3(c,d)]. This per- 
turbed R292 conformation is also observed in a calcium 
deficient crystal structure of N9 neuraminidase. In the 
AMBER FF99SB ion-free simulations, interaction with 
oseltamivir and R292 is weakened (Table I), but this con- 
formational change is not observed. 

Binding free energies 

To quantify the impact of the aforementioned changes 
in dynamics of the Nl-oseltamivir complex in absence of 
calcium, parallel free energy calculations were performed: 
free energy perturbation (FEP) with the AMBER FF99SB 
force field, post-processed with the MBAR method, and 
IT-TI with the GROMOS96 force field. 

For comparison to the ion-bound FEP results (Table 
III) and to investigate the contribution of Y347 to the 
binding free energy, two sets of free energy calculations 
for the ion-free tetramer were carried out. The first cal- 
culation was initialized from a calcium-free structure 
selected after 5 ns of equilibration. During this equilibra- 
tion period, Y347 had not yet flipped out of the pocket 
for three of the four monomers (see Supporting Informa- 
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a 



N294 




Y347 



I222 



oseltamivir 



N294 



R224 



I222 



Y347 



E276 

Is- 

C177 



oseltamivir 



Figure 3 

Binding pose comparison for ion-bound and ion-free simulations. Oseltamivir is shown in dark gray for all panels, with residues from the top 
cluster of the GROMOS96 ion-bound simulation (pink) in panels (a) and (c), from the cluster four of the GROMOS96 ion-free simulation 
(yellow) in panel (b), and from the cluster two of the AMBER FF99SB ion-free simulation (cyan) in panel (d). The disturbed hydrophobic 
subpocket structure is shown with surface representation in (b) and the retreated R292 conformation in (d). 



tion Fig. 2) and monitoring of the torsion throughout 
this calculation indicated 17% occurrence of the "out" 
Y347 state, compared to 6% in the calcium-bound simu- 
lation (black solid vs. black dashed in Fig. 4). Here, the 
free energy difference for the ion-bound and ion-free 
complexes is very small (Table III), which suggests that 
the AMBER FF99SB force field confers stability to other 
binding site residues and preserves an effective binding 
pose with oseltamivir. However, when cluster two (Table 
II) of the AMBER FF99SB ion-free MD simulations was 
selected as starting structure for each monomer in a sec- 
ond tetramer calculation, simulation time spent in the 
Y347 "out" conformation increased to 76% (green dashed 
line, Fig. 4). This increased destabilization of the active 
site framework gives more impact on the free energy of 
binding, with a +2.7 kcal/mol reduced favorability in 
AAGbind- The change in results underscores both the 
influence of starting structure on sampling in free energy 
calculations, as well as the difficulty in achieving complete 
sampling using traditional MD methods. 

The TI results are an average of six monomer calcula- 
tions, initiated by structures that capture diverse confor- 



mations of Loop 150, as well as varied conformations of 
Y347. For the ion-free complexes, three of the starting 
structures include a Y347 "out" conformation, and the 
population of this conformation increases from 13% to 
20% in the ion-free calculations (Fig. 4). This population 
shift is not as significant as seen for the GROMOS96 
MD (Fig. 2), due to the shorter timescale of the sampling 
in each TI X window and fewer independent trajectories 
(6 for IT-TI, 10 for the MD). Still, the calcium-free 



Table III 

Binding Free Energies (kcal/mol) 





lon- 


bond 


Ion-free 






AGbind 


Tffi 


AG bind 


a a 
\ N 


AMBER FF99SB with Y347 in 


-14.5 


1.5 


-14.3 


0.7 


AMBER FF99SB with Y347 out 






-118 


1.3 


GR0M0S96 


-15.0 


2.1 


-10.4 


1.2 


Experiment* 1 


-13.7 


0.2 


n/a 




a JV = 4 for AMBER calculations 


and N — 6 


for GROMOS 


calculations 


(see 



Methods). 

b Derived from iCj value by Kati et al.^ 
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Figure 4 

Y347 torsion changes during free energy calculations for the AMBER 
FF99SB ion-bound calculation (black, solid lines), AMBER FF99SB ion- 
free calculation with Y347 "in" starting structure (black, dashed), and 
with Y347 "out" starting structure (green, dashed); Also, GROMOS96 
ion-bound simulations (red, solid lines) and GROMOS96 ion-free 
simulations (red, dashed). 



AGbind is +4.6 kcal/mol less favorable than the ion- 
bound free energy value (Table III). Structural evidence 
for the larger impact of the calcium presence on the 
GROMOS96 free energy compared to the AMBER 
FF99SB results can be found in the more pronounced 
disruptions of the active site framework, a more signifi- 
cantly perturbed binding pose, and overall more flexibil- 
ity observed in the GROMOS96 MD. However, the dis- 
crepancy between the results for the two force fields is 
within the reported, propagated error. 

The free energy results, listed in Table III, are a close 
match to each other and reflect a very high, nanomolar 
binding affinity. From a survey of neuraminidase activity 
assays in the literature^ - " j- ne average IC 50 value for 
oseltamivir binding to Nl is ^0.2 nM. In Table III, we 
compare the calculated results with a free energy value 
derived from the Nl -oseltamivir Kj in a kinetic study by 
Kati et al. Both free energy results are close to this 
value, within the reported error. Because all binding 
assays in the literature include calcium salts in the activ- 
ity buffer, we lack a direct experimental comparison for 
the ion-free results, but the reduced favorability, seen in 
the GROMOS free energy results and the AMBER results 
for predominantly "out" Y347 conformations (Table III), 
corresponds to an increase in K; by ~2 orders of magni- 
tude (from AGb ind = +RT\nKi). 

The free energy calculations, considered with the struc- 
tural changes observed in the MD simulations, implicate 
calcium as an important factor in maintaining an effec- 
tive Nl active site framework, and, in particular, through 
stabilization of Y347. The "in" conformation of this resi- 
due seems to "clamp" the ligand into a favorable binding 
pose, which includes consistent interactions with R292 
and the rest of the R triad, as well as orientation of osel- 



tamivir such that its branched aliphatic chain is directed 
to a nearby hydrophobic subpocket. The sensitivity of 
neuraminidase binding affinity to changes in this region 
of the active site is also evidenced by the oseltamivir- re- 
sistant mutations of structurally neighboring residues: 
R292K (N2), N294S (Nl), and H274Y (Nl). 4 Y347 has 
been suggested ^ to prevent oseltamivir resistance to the 
R292K mutation in most Nl strains by providing extra 
stability to the binding pose, as is observed in our simu- 
lations and supported by free energy calculations. How- 
ever, Y347 is mutated to an asparagine in the Nl protein 
of both the 1918 epidemic influenza strain and one of 
examined 2009 "swine flu" strains,^ which could allow 
susceptibility for resistance mutations. 



CONCLUSION 

Consistent structural changes observed in simulations 
with two different force fields give reinforced evidence 
for the structural role of calcium in stabilization of the 
Nl active site. Free energy results calculated using differ- 
ent methods and force fields agree closely with each 
other and with experiments to indicate a contribution of 
3-5 kcal/mol to the binding free energy, or two orders of 
magnitude in _K";. The methods in this study also used 
multiple, independent simulations, initiated from diverse 
starting structures, in an effort to enhance sampling of 
the MD and obtain more accurate and reliable free 
energy results. Even on timescales as short as 1 ns 
(GROMOS96 TI), absence of the calcium causes disrup- 
tion of the active site, particularly in the flipping "out" 
of Y347, and perturbation of the ligand binding pose. In- 
tegrity of the active site is essential for structural analysis 
and modeling of the Nl-drug interface, often used in 
rational drug design. We emphasize the importance of 
the calcium for maintenance of this integrity and urge 
that it should not be excluded in Nl active site studies, 
particularly for structures used in docking and other 
computational methods. For the challenges of emerging 
resistance mutations in Nl, which have so far occurred 
in the region of the active site near the bound calcium, 
noting the importance of framework residues in this area 
could also guide design of new, robust inhibitors. 
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